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Abstract 

We present a new algorithm for isothermal-isobaric molecular-dynamics simulation. The method uses an 
extended Hamiltonian with an Andersen piston combined with the Nose-Poincare thermostat, recently developed 
by Bond, Lcimkuhler and Laird [J. Comp. Phys., 151, xxxx (1999)]. This Nose-Poincare-Andersen (NPA) 
formulation has advantages over the Nose-Hoover- Andersen approach in that the NPA is Hamiltonian and can 
take advantage of symplectic integration schemes, which lead to enhanced stability for long-time simulations. The 
equations of motion are integrated using a Generalized Leapfrog Algorithm and the method is easy to implement, 
symplectic, explicit and time reversible. To demonstrate the stability of the method we show results for test 
simulations using a model for aluminum. 

, 1 Introduction 

> ■ . . 

q>^ Traditionally, molecular-dynamics simulations are performed using constant particle number N, volume V and 
energy E. However, these are not usually the conditions under which experiments are done and there has been much 
<0 attention to the development of simulation methods designed to sample from other, experimentally more relevant 
C") ensembles, such as constant temperature (canonical) and/or constant pressure^, |j. Some of the most popular and 
useful of these are those based on so-called "extended" Hamiltonians, i.e., Hamiltonians in which extra degrees of 
freedom have been added to the system in order to ensure that the trajectory samples from the statistical distribution 
*"^» corresponding to the desired thermodynamic conditions. 

For a constant pressure system, for example, Andersen^] introduced the volume V, along with its corresponding 
conjugate momentum ny, as extra variables. The new variables are coupled to the system in such a way as to 
^jguarantee that the trajectory (if ergodic) samples from an isobaric statistical distribution. Similarly, to generate a 
,£h constant temperature distribution, Nosey] introduced a new mechanical variable s (with conjugate momentum 7r s ) 
CUhat couples into the system through the particle momenta and acts to effectively rescale time in such a way as 
to guarantee canonically distributed configurations. These two extensions can be combined to give a Hamiltonian 
whose trajectories can be shown to sample from an isothermal-isobaric ensemble. M 
This combined Nose- Andersen (NA) Hamiltonian is given by 
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where pi is the conjugate momentum to the scaled position qi = V - 1 ' i Vi, P ex t is the external pressure and g is 
given by Nf + 1 where Nt is number of degrees of freedom of the original system. The quantities Qy and Q a are 
the masses of the Andersen "piston" and the Nose thermostat variable, respectively. 

The equations of motion for this system are 

Pi = -V^ViUiV^q) (2a) 
^ Pi (2b) 
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TTV = P-Pext (2c) 

V" = 7r,/Qy (2d) 

tt s = ^ 2 / 3 S - 3 V^-^ (2e) 

a = Ks/Qs, (2f) 



where the instantaneous pressure V is given by 
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3^ ^ 2m l V 2 / 3 s 2 3V ^ dq. 

There are two major drawbacks to this approach: First, because of the time rescaling, the time variable in Nose 
dynamics is not "real" time, so any discretized trajectory generated by numerically integrating the Nose equations 
of motion must be transformed back into "real" time, leading to the configurations that are spaced at unequal 
"real" -time intervals. This is inconvenient for the construction of equilibrium averages, especially of dynamical 
quantities. Second, the Hamiltonian is not separa6/e||] (that is, the kinetic and potential terms in the Hamiltonian 
are not functions only of momenta and position variables, respectively), making standard Verlet/leapfrog approaches 
inapplicable. 

By a change of variables and a time rescaling of the equations of motion, Hoover || derived new equations of 
motion that generate the same trajectories (for the exact solution) as the original Nose Hamiltonian, but in real time. 
This Nose- Hoover dynamics has become a standard method in molecular simulation. However, the change of variables 
that links the Nose Hamiltonian to the Nose-Hoover equations of motion is a non-canonical transformation - the 
total energy function of the system is still conserved, but it is no longer a Hamiltonian, since the equations of motion 
cannot be derived from it. Although a variety of very good time-reversible methods have been put forward, |To|, O] 
the lack of Hamiltonian structure precludes the use of symplectic integration schemes, which have been shown to 
have superior stability over non-symplectic methods 

2 The Nose-Poincare- Andersen (NPA) Hamiltonian 

Recently, Bond, Leimkuhler and LairdjlJ] have developed a new formulation of Nose constant-temperature dynamics 
in which a Poincare time transformation is applied directly to the Nose Hamiltonian, instead of applying a time 
transformation to the equations of motion as in Nose-Hoover. The result of this is a method that runs in "real" time, 
but is also Hamiltonian in structure. In this work we combine this new thermostat with the Andersen method for 
constant pressure to give an algorithm for isothermal-isobaric molecular dynamics. For a system with an Andersen 
piston, the new Nose-Poincare-Andersen (NPA) Hamiltonian is given by 

Hnpa = [H NA - n NA {t = 0)}s , (4) 

where Hn a is given in Eq. [I]. As discussed in Ref. 12, the above form of the Hamiltonian (a specific case of a 
Poincare time transformation) will generate the same trajectories as the original Nose- Andersen Hamiltonian, except 
with time rescaled by s (which puts the trajectories back into real time), The resulting equations of motion (except 
for 7r s ) for this constant pressure and temperature Nose-Poincare Hamiltonian are the as those given above for the 
Nose-Andersen system (Eqs. pa}|2f|), except that the right-hand side is multiplied by the thermostat variable s. For 
7r s we have 

n s = -s^-AH = V- 2 / 3 Y ^L- gkT -AH (5) 

OS ^— ' TfliS 

where AH = Una - Hna^ = 0). 



3 Integrating the NPA Equations of motion 

The NPA Hamiltonian is nonseperable since the kinetic energy contains the extended "position" variables s and 
V. The equations of motion for a general time- independent, non-seperable Hamiltonian can be written (for general 
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positions Q and conjugate momenta P) 



Q = G(P,Q) 

P = F(P,Q), (6) 

where G(P, Q) — and F(P, Q) = — (For a seperable Hamiltonian G is only a function of P and F is only a 
function of Q.) For such a nonseperable system, standard symplectic splitting methods, such as the Verlet/leapfrog 
algorithm, are not directly applicable. However, symplectic methods specifically for nonseperable systems have been 
developed ||. One simple example that is second-order and time-reversible is the Generalized Leapfrog Algorithm 
(GLA) 

Pn+l/2 = Pn + hF{P n+1/2 ,Q n )/2 
Qn-\-l — Qn 

+ h[G(P n+1/2 ,Q n ) + G(P n+ 

1/2 j Qn+1 )]/2 

Pn+l = Pn+l/2 + hF(P n+1/2 ,Q n+1 )/2, (7) 

where h is the time step and P n and Q n are the approximations to P(t) and Q(t) at t = t n = nh. (This method can 
be obtained as the concatenation of the Symplectic Euler method, 

Pn+l = Pn + hF(P n+ i,Q n ) 
Qn+1 — Qn 

+ hG(P n+1 ,Q n ) , (8) 

with its adjomtQ. The concatenation of a integrator with its adjoint guarantees a time-reversible method.) This 
method is a simple example of a class of symplectic integrators for nonseperable Hamiltonians[jl3[ [l4], |l5|, |l6| |. 
Applying the GLA to the NPA equations of motion gives 

Pi,„+l/2 = P»,r l -^„y„ 1 / 3 V l C/(^ 1 / 3 q) (9a) 

7T v,n+l/2 — 7: vn + —s n [P(q n ,p n+ i/2,V n ,s n )—P ext \ (9b) 

h I A Pjn+1/2 , \ 

\i=l W 8 Vn s„ / 

— ■^■^W(qn: Pn+l/2> V n ,T7 v n+ i/2, S n , T^ s ,n+l/2) (9°) 
S„+l = S n + -(s„ + S n +l) s ' n + 1/2 (9d) 

Vn+1 = K + -(s n + s„ + i) -pr— — (9e) 

, ^ / 1 , 1 \ Pi,n+l/2 , Q ~ 

= g -+2 7^ + (9f) 



v S„lV S„+iV„ +1 
7Ts,„ + l = ^+1/2+2 ^ /3 - 9 k B T\ 

- 2^^( Ct ™+ 1 ' Pn+l/2) 7r «,n+l/2, Sn+l, 7T S ,n+l/2) (9g) 

Kv,n+1 = ^v,n+l/2 + -^S n+ i[P(c[ n+ i,p n+1 /2,V n+ i,S n+ i) - P ext ] (9h) 

1/3 1/3 

Pi,n+i = iV+i/2 + 2 Sn + lV n+i ViU(V n+1 q. n +i) (9i) 



As in the case of the constant volume Nose-Poincare algorithm, the GLA for the NPA is explicit - this is not 
necessarily the case for a general nonseperable Hamiltonian. Note that Eq. ^c] requires the solution of a scalar 
quadratic equation for tt s >n +i/2- Details of how to solve this equation without involving subtractive cancellation can 
be found in Ref. 12. 
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4 Test Simulation Results and Summary 



In order to evaluate this method, simulations were performed using an embedded atom potential for aluminum |L7j. 
We report test simulations on a system of 256 particles with periodic boundary conditions for an aluminum melt at 
T = 1000-RT and P = 0. For this model mass is measured in amu, distance in A and energy in eV, the natural time 
unit of the simulation is then 10.181fs, that is, a simulation time step of 0.1 corresponds to an acutal time step of 
1.0181fs. 

First the stability of the method was tested. Fig. 0(a) shows the value of the NPA Hamiltonian (a conserved 
quantity) as a function of time in a long run. The trajectory shown here was begun after initial equilibration at 
1000K for 2xl0 6 time steps (2.03ns). In this simulation, the values of Q v and Q s (in reduced units) were 10~ 4 
and 2.5, respectively. The stability of the method is excellent, giving no noticeable drift in Hnpa over ^ e course 
of a long trajectory. The pressure and temperature trajectories for this run are also shown in Figs. 0(b) and 0(c), 
respectively. 
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Figure 1: (a) The value of the NPA Hamiltonian as a function of time for a long run (over 2ns) on a 256 particle aluminum system. 
The starting configuration had been equilibrated at T = 1000K and P = 0. The values of Q v and Q s (in reduced units) are 0.4 and 
10 — 4 , respectively. The pressure and temperature trajectories for this run are given in (b) and (c), respectively. 



In Fig. 2(a)-(c), we show the ability of the method to regulate the pressure, temperature and density, respectively, 
for three sets of extended variable masses. (The values for Q s were larger here than that used in Figure p] because 
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small values of that variable lead to instabilities when the initial temperature is far from the target temperature.) 
The system was initialized to an fee lattice of initial density 0.0602lA~ 3 with the individual velocity components 
chosen from a Maxwell-Boltzmann distribution at 100K. The simulation was then run with the NPA with T — 1000A" 
and P = 0. In all cases, the instantaneous pressure, temperature and density evolve quickly and stabilize about their 
desired values. 
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Figure 2: Pressure (a), temperature (b) and density (c)trajectories for the 256 particle aluminum system using the NPA 
algorithm with T = 1000if and P = starting from an initial configuration in an fee lattice with density p — 0.06021A -3 and 
initial velocities chosen from a Maxwell-Boltzmann distribution at 100K. 



The GLA has a global error that is second-order in the time step. To demonstrate that this also is true in our 
results, a series of simulations were performed using various values for the time step. The system was initialized in 
an identical manner to that described in the last paragraph and then run for a total time of 2.0362 ps. Figure S 
shows, for several combinations of extended variable masses, a log-log plot of the energy error, as estimated by the 
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standard deviation of Ti-^pj^ (Eq. versus the time step. One notes that here smaller values of Q v lead to smaller 
fluctuations in Hnpa.. 





Figure 3: Log- log plot of the energy error STLnpa versus the time step for a variety of piston and thermostat masses, 
order of the error is given by the slope of lines and is labeled as m in the legend. 
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Figure 4: The normalized velocity autocorrelation function, C(t), for the embedded-atom model for aluminum calculated for 
a 256 particle system at T=1000K and P=0. The solid line is calculated using the NPA algorithm described herein. The 
dotted line is a constant NVE simulation under similar conditions. 

Finally, to demonstrate that the method yields relevant dynamical quantities, the normalized velocity autocor- 
relation function, C(t) = (v(i) • v(0))/(v(0) • v(0)), was calculated using our constant NPT algorithm (with Q v and 
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Q s as in Fig. 1) and compared to the same quantity calculated using standard constant NVE molecular dynamics 
(with a velocity- Ver let integrator ]L8|). The NVE simulations were run at an energy and density corresponding to 
the average energy and density for the constant NPT simulations. This comparison is shown in Figure [Q. Both 
systems were first equilibrated at 1000K for 200,000 steps (203.6 ps) and run for 20,000 steps (20.36 ps) to collect 
averages. C(t) for the Nose-Poincare- Andersen method for constant NPT molecular dynamics is indistinguishable in 
this figure from that of the NVE simulation. 
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